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Abstract 

For asymptotically periodic systems, a powerful (phase) reduction of the dynamics is obtained by 
computing the so-called isochrons, i.e. the sets of points that converge toward the same trajectory 
on the limit cycle. Motivated by the analysis of excitable systems, a similar reduction has been 
attempted for non-periodic systems admitting a stable fixed point. In this case, the isochrons 
can still be defined but they do not capture the asymptotic behavior of the trajectories. Instead, 
the sets of interest — that we call "isostables" — are defined in literature as the sets of points that 
converge toward the same trajectory on the stable slow manifold of the fixed point. However, it 
turns out that this definition of the isostables holds only for systems with slow-fast dynamics. Also, 
efficient methods for computing the isostables are missing. 

The present paper provides a general framework for the definition and the computation of the 
isostables of stable fixed points, which is based on the spectral properties of the so-called Koopman 
operator. More precisely, the isostables are defined as the level sets of a particular eigenfunction of 
the Koopman operator. Through this approach, the isostables are unique and well-defined objects 
related to the asymptotic properties of the system. Also, the framework reveals that the isostables 
and the isochrons are two different but complementary notions which define a set of action-angle 
coordinates for the dynamics. In addition, an efficient algorithm for computing the isostables is 
obtained, which relies on the evaluation of Laplace averages along the trajectories. The method 
is illustrated with the excitable FitzHugh-Nagumo model and with the Lorenz model. Finally, 
we discuss how these methods based on the Koopman operator framework relate to the global 
linearization of the system and to the derivation of special Lyapunov functions. 

Keywords: Nonlinear dynamics, isochrons, excitable systems, Koopman operator, action-angle coordinates, 
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I. INTRODUCTION 



Among the abundant literature on networks of coupled systems, a vast majority of studies 
focus on asymptotically periodic systems (i.e. coupled oscillators) while only a few consider 
coupled systems characterized by a stable fixed point. This is particularly surprising since the 
atter can exhibit excitable regimes that are relevant in many situations (e.g. neuroscience 



14J). One reason for this disproportion is probably related to phase reduction methods. 



For asymptotically periodic systems, powerful phase reduction methods turn the (complex, 
high-dimensional) system into a phase oscillator evolving on the circle, making the analysis 



3JJ. In contrast, in 



of complex networks more amenable to mathematical analysis 

an 

the case of systems admitting a stable fixed point, the development of equivalent reduction 
methods is more recent and a general framework is still in its infancy. 

The goal of reduction methods is to assign the same value to a (codimension-1) set of 
initial conditions that are characterized by the same asymptotic behavior, in turn designing 
a coordinate on the state space. In the case of asymptotically periodic systems, these sets of 
identical (phase) yalue are the so-called isochrons, which approach the same trajectory on 



the limit cycle 



321 ] . This concept has been recently extended to heteroclinic cycles 29( . For 



systems admitting a stable focus, the isochrons (or isochronous sections) can still be defined 

nn 

as the sets of points that are invariant under a particular return map [9|, [28[ . This notion is of 
particular interest in the case of weak foci (i.e. with purely imaginary eigenvalues) and non- 
smooth vector fields, where the existence of isochrons is a non-trivial problem related to the 
stability of the fixed point. However, the isochrons provide in this case no information on the 
asymptotic convergence of the trajectories toward the fixed point and are not useful for the 
system reduction. (Note also that they do not exist for fixed points with real eigenvalues.) 
Therefore, the isochrons must be complemented by another family of sets: the so-called 
isostables. 

Excitable systems are characterized by slow-fast dynamics with a stable fixed point and, 
in the plane, they admit a particular trajectory — the transient attractor or slow manifold — 
that temporarily attracts all the trajectories as they approach the fixed point. In this case, 
the isostables are naturally defined as the sets of points that converge to the same trajectory 



on the transient attractor 



24l | . (Note that these sets are called "isochrons" in 24j, but we 



feel that the proper sense is "isostables" instead, in order to avoid the confusion with the 
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isochrons of foci studied 



in 0, 



|28j.) For non-planar systems possessing a multi-dimensional 



slow manifold or center manifold, a (more rigorous) framework was previously developed 



in 



26]. In that work, the sets of interest (called "projection manifolds" in 26]) are 



closely related to the notion of isostable and correspond to the invariant fibers of the (slow 
or center) manifold, i.e. the sets of initial conditions characterized by the same long-term 
behavior on that manifold. Through the reduction obtained with the isostables, excitable 
systems have been studl ed b various con te * ts (sensitive to period* pulses £ Q Q, 
network synchronization 17j, etc.). 

Since the isostables provide a characterization of the system dynamics around the fixed 
point, their computation is also desirable for systems which do not contain multiple time 
scales (i.e. with no slow or center manifold). For instance, the computation of the isosta- 
bles can be useful to achieve an optimal control that minimizes the time of convergence 
toward a steady state or to investigate the delay of convergence to a stable equilibrium 
in decision-making models 3(|. But in these cases, a more general framework is required, 
which defines the isostables as particular (and unique) codimension-1 sets capturing the 
asymptotic behavior of the system. In addition, the computation of the isostables through 



backward integration 



24[ | or normal form of the dynamics 



4[ is limited to a neighborhood 



of the slow manifold. In this context, an efficient method for computing the isostables in 
the entire basin of attraction is also missing. 

In this paper, we propose a general framework for the reduction of systems admitting a 



stable fixed point, which is not limited to excitable systems with slow-fast dynamics. 



'his 



approach is based on the spectral properties of the so-called Koopman operator [19|, [21]. 
More precisely, we propose a general and unique definition of the isostables in terms of a 
particular eigenfunction of the Koopman operator. In addition, the framework yields an 
efficient method to compute the isostables in the whole basin of attraction. This method 
relies on the estimation of Laplace averages al ong the trajectories and can be seen as an 



extension of the approach recently developed in [18] to compute the isochrons of limit cycles. 

Viewed through the Koopman operator framework, the isostables and the isochrons ap- 
pear to be two different but complementary concepts. On the one hand, they are different 
since they are related to the absolute value and to the argument, respectively, of the eigen- 
function of the Koopman operator. On the other hand, they are complementary in the sense 
that they define a set of action-angle coordinates for the system dynamics. This action-angle 
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representation is related to important properties of the isotables, such as the global lineariza- 
tion of the dynamics and the derivation of special Lyapunov functions, that we discuss in 
the paper. 

The paper is organized as follows. In Section [Til we introduce the concept of isostable in 
the context of the Koopman operator framework, both for linear and nonlinear systems. We 
also propose a rigorous definition of the isostables and discuss their main properties. The 
relation between the isostables and the Laplace averages is developed in Section IIHI This 
provides an efficient algorithm for the computation of the isostables which is illustrated in 
Section IIVI for the excitable FitzHugh-Nagumo model and the Lorenz model. Finally, the 
related concepts of action-angle representation, global linearization, and Lyapunov function 
are discussed in Section [V] Section I VI I gives some concluding remarks. 

II. ISOSTABLES AND KOOPMAN OPERATOR 

The isostables of an asymptotically stable fixed point x* are the sets of points that share 
the same asymptotic convergence toward the fixed point. More precisely, trajectories with 
an initial condition on an isostable X ro simultaneously intersect the successive isostables X Tn 
after a time interval r n — To, thereby approaching the fixed point synchronously (Figure 
[TJ). The isostables partition the basin of attraction of the fixed point and define a new 
coordinate r that satisfies f = 1 along the trajectories. Or equivalently, they define a 
coordinate r = exp(Ar) with the linear dynamics r = \r. This new coordinate can be used 
in a context of model reduction. 

At this point, it is important to remark that this (intuitive) definition of isostable is not 
complete. Indeed, there exist an infinity of families of sets that satisfy the above- described 
property. In this section, we will show that it is possible to give a rigorous definition of 
the isostables, which are unique and well-defined objects. To do so, we first consider the 
particular case of linear systems. Then, we extend the concept to nonlinear systems, using 
the Koopman operator framework. 
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Figure 1: Trajectories starting from the same isostable I TQ are characterized by the same con- 
vergence toward the fixed point. They simultaneously intersect the successive isostables l Tn and 
approach the fixed point synchronously. 

A. Linear systems 

Consider the stable linear system 

x = Ax, x G R n , (1) 

and assume that each eigenvalue Xj = <x, + iujj of the matrix A is of multiplicity 1, has 
a strictly negative real part Uj < 0, and corresponds to the right eigenvector Vj (which is 
normalized, that is, ||vj|| = 1). By convention, we sort the eigenvalues so that Ai is the 
eigenvalue related to the "slowest" direction, that is 

<Tj < o-i < , j = 2, . . . , n . (2) 

The flow induced by (CQ) is the continuous-time map : K x M. n i— > W 1 , that is, <fr(t, x) is the 
solution of (p]) with the initial condition x. For linear systems, the flow is given by 

n 

0(t,x) = ^ Sj (x) Vj e^, (3) 
i=i 

where Sy(x) are the coordinates of the vector x in the basis (v 1; . . . , v n ). The function Sj(x) 
can be computed as the inner product Sj(x) = (x, Vj), with Vj the eigenvectors of the adjoint 
A*, associated with the eigenvalues = crj —iu>j and normalized so that (v^, v\,-) = 1. (Note 
that Sj(x) is an eigenfunction of the so-called Koopman operator; see Section Til Bl ) 

Next, we show that the isostables of linear systems are simply defined as the level sets of 
|si(x)| = |(x, Vi)j. We consider separately the cases Ai real and Ai complex. 
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1. Real eigenvalue \\ 



When the eigenvalue Ai = <j\ is real, the trajectories induced by the flow ([2D asymp- 
totically approach the fixed point along the slowest direction vi (since the eigenvalues 
are sorted according to ((2D). Then, the trajectories characterized by the same coefficient 
|si(x)| = exp(o"ir(x)) exhibit the same asymptotic convergence toward the fixed point: 

n 

<j>{t, x) = v±e ffl( * +r(x)) + J2 Sj(x) v, exp(\jt) « v ^e <Tl(t+T(x)) as t -)• oo , (4) 

i=2 

where the notation vf implies that either the vector vi or — vi must be considered. The 
initial conditions x of these trajectories therefore belong to the same isostable 



x G 



x = e CTir v± + £ a, v j , Wa, G 
i=2 



(5) 



which is obtained by considering t — in (jlj. In this case, the isostables are the (n — 1)- 
dimensional hyperplanes parallel to v 3 - for all j '• > 2 (or equivalently, perpendicular to vi) 
(Figure [2D- 





(a) 



(b) 



Figure 2: (a) The isostables of linear systems with a real eigenvalue Ai are the hyperplanes spanned 
by the eigenvectors v,-, with j > 2. The particular isostable 1^ contains the fixed point, (b) For 
two-dimensional systems (or in the plane vi — V2), the isostables correspond to pairs of parallel 
lines. 
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2. Complex eigenvalue Ai 



A system having a complex eigenvalue Ai can be transformed through the use of action- 
angle coordinates. Then, the isostables are obtained from the isostables fl^D of the subsystem 
which is related to the action coordinates and which is only characterized by real eigenvalues 
Oy Consider a linear coordinate transformation that expresses the dynamics (pQ) in the 
(spectral) basis given by the vectors Vj (for Aj real) and 3ft{vj}, — 3{vj} (for Aj = Aj +1 
complex). (Note that 9?{vj} and ^{vj} are not parallel since the two eigenvectors Vj 
and Vj+i are independent.) This is performed by diagonalizing A and by using the linear 
transformation 



/ 

V 



1 1 

-i i 



in each subspace spanned by a pair of complex eigenvectors (vj,Vj + i). The dynamics become 



Vj 
( 



Vj 



\ ( 



Xj real , 



-OJn 



( 



3 ^3 



\ w i j J 



\ 



Vj 



A; 



Aj +1 complex . 



with the initial conditions Hj(0) 
(23fJ{ Sj (x )},23{ Sj (x )}) (for A, = X c j+1 
Xj G R) and the polar coordinates (yj, Uj+i) 
obtain the canonical equations 



Sj(xo) (for Aj G 



and (yj(0),y j+1 (0)) 



I). Then, using the variables rj = yj (for 
i r j cos(^), rj sin(^)) (for \ = X C J+1 £ R), we 



°3 T 3 > 



U)j. 



(6) 
(7) 



The initial conditions are given by r,(0) = Sj(x ) (for Aj G R) and (^j(O), 0j(0)) = 
(2| S j(x )|,Z S j(x )) (for A, = X< +1 $ R). 

According to (EJ)— C[3) , the variables Vj and 9j can be interpreted as the action-angle co- 
ordinates of the system (see [lj]) and the convergence toward the fixed point is captured by 
the (action) variables rj. Therefore, the isostables of ([I]) correspond to the isostables of the 
affine system with the real eigenvalues aj. Given the results of Section Hi A 1\ it follows 
that the isostables are characterized by a constant value |ri|, that is, they are the level sets 
of |si(x)|. Denoting r\ = 2|si(x)| = exp^i^x)) and using an expression similar to ([3]), we 
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obtain (in the variables y. 

?r = jy e 

where are the unit vectors of M n , or equivalently (in the variables x 



y = (cos(0)ei + sin(0)e 2 )e <Jir + £ a, e, , Va, £ R, 6 [0, 2tt) 

i=3 



x G 



x = (cos(#)a + sin(^)b)e <T1T + ^ v,- , Voij G K, V6* G [0, 2tt) ) , (8) 

i=3 



with a = 5i{vi} and b = — 3{vi}. In this case, the isostables are the (n — l)-dimensional 
cylindrical hypersurfaces parallel to v, for all j > 3. The intersection of an isostable with 
the 2- dimensional plane spanned by (a, b) (i.e., the base of the cylinder) is an ellipse (Figure 
Indeed, a linear transformation turns the circle in the variables yj into an ellipse in the 
variables Xj. 

The trajectories starting from the same isostable converge to the fixed point along a spiral 
characterized by the vectors (a, b), according to 

x) » (a cos(wit + 0(x)) + b ssnfat + 0(x))) e ai(t+T{x)) as t ->• oo , 



with exp(o"ir(x)) = 2|si(x)| and 9(x) = Zsi(x). Note that the phase — or angle coordinate 



9 is related to the notion of isochron (see e.g. 




Q,|28| and Section Wj- 

/ 9 = 7T/2 



(a) 




0(t,x) 



Figure 3: (a) The isostables of linear systems with a complex eigenvalue Ai are cylindrical hyper- 
surfaces spanned by Vj for all j > 3. (b) For two-dimensional systems (or in the plane a — b), the 
isostables correspond to ellipses. 



The expressions and (|SJ) provide a unique definition of the isostables in the case of 
linear systems, when Ai is real and when Ai is complex, respectively. Since vf = Vi exp(z#) 
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with 6 = {0,7r} and cos(#)a + sin(0)b = 3l{vi exp(i8)}, these two definitions can be sum- 
marized in a single definition. 



Definition 1 (Isostables of linear systems). For the system ([I]), the isostable I T associated 
with the time r is the (n — lVdimensional manifold 



1 T = { x G £(x*; 



x = |vi e ie } e CTlT + ^ ay v,- , Va, el,V^9 



3=3 



with 6 = {0, tt} and j = 2 if \ e R, and 6 = [0, 2vr) and j = 3 if A x £ R. 

B. Nonlinear systems 

Now, we consider a nonlinear system 

x = F(x) , x6R" (9) 

where F is an analytic vector field, which admits a stable fixed point x* with a basin of 
attraction £>(x*) C M. n . In addition, we assume that the Jacobian matrix J computed at 
x* has n distinct (nonresonant) eigenvalues Xj = a,j + iojj characterized by strictly negative 
real parts Oj < and sorted according to (j2J). (For unstable fixed points or for multiple 
eigenvalues, see Remark [T] and Remark [2J respectively.) 

The isostables of linear systems have been defined as the level sets of the coefficient si(x) 
that appears in the expression of the flow ([2]). For nonlinear systems, an expression of the 



flow similar to ([3D can be obtained through the framework of Koopman operator 19|,|2l|. The 
Koopman semigroup of operators U l describes the evolution of a (vector-valued) observable 
f : W l H> C m along the trajectories of the system and is rigorously defined as the composition 
t/*f(x) = f o0(t,x). In particular, for an analytic observable, the spectral decomposition of 



the operator yields 



20] 



C7*f(x) = £ ^(x)---^"(x)v fcl ... fc?i e( fclAl+ - +fc " A ") i . (10) 
{fci,...,fc„}eN™ 

A detailed derivation of the decomposition in the case of a stable fixed point is given in Ap- 
pendix [A] The functions Sj(x), j — 1, . . . , n, are the smooth eigenf unctions of U l associated 
with the eigenvalues Xj, i.e. 

f/S-W = ^(^,x)) = Sj (x)e^, (11) 
10 



and the vectors v fcl ... fcn are the so-called Koopman modes [27(, i.e. the projections of the 
observable f onto s 1 1 (x) • • ■ s^ n (x). For the particular observable f(x) = x, (ll(J|) corresponds 
to the expression of the flow and can be rewritten as 



(f)(t, x) = C/*x = x* + s i( x ) v i e 



A,/ 



E 



W v fcl ... fcn e 



(feiAiH hfc„A„)t 



{fei,...,A: n }€N™ 
fel+-+fc«>l 



(12) 



The first part of the expansion is similar to the linear flow ([3]). The eigenvalues Xj and the 
Koopman modes Vj are the eigenvalues and eigenvectors of J, respectively. Although the 
eigenfunctions s, (x) are not computed as the inner products (x, Vj) as in the linear case, 
they can be interpreted as the inner products (z,Vj), where z is the initial condition of a 
virtual trajectory evolving according to the linearized dynamics z = Jz and characterized 
by the same asymptotic evolution as </>(£, x) 3]. The other terms in (I12p do not appear in 
the expression of the linear flow and account for the transient behavior of the trajectories 
owing to the nonlinearity of the dynamics. 

The isostables can be rigorously defined as the level sets of the absolute value of the 
eigenfunction |si(x)|. Indeed, the asymptotic evolution of the flow f[T2l is dominated by 
the first mode associated to Ax- Then, a same argument as in Section Til Al shows that the 
points x characterized by the same value |si(x)| are the initial conditions of trajectories that 
converge synchronously to the fixed point, with the evolution 



0(t,x) 



x* + v 1 t e Al( * +T(x)) , e CTir(x) = |si( 

x * _|_ Sft j Vl e «(wi*+0(x)) j e CTi(t+r(x)) ^ 



x)|, Ai G 

e aiT(x) = 2|si(x)| , 0(x) = Zsi(x) X 1 i 



(13) 



We are now in position to propose a general definition for the isostables of a fixed point, 



which is valid both for linear and non 
definition of isochrons for limit cycles 



inear systems and which is reminiscent of the usual 
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22]. 



Definition 2 (Isostables). For the system (JUJ), the isostable X T of the fixed point x*, asso- 
ciated with the time r, is the (n — lVdimensional manifold 



1 T = x G #(x*; 



36eQs.t. lime - ' 71 * 

t— »oo 



0(t,x) -x* - 3?{v 1 e i(a;it+e) }e <Tl(<+T) || = o| , 



with 6 = {0, ?r} if Ai e R and 6 = [0, 2vr) if A x i R. 
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The reader will easily verify that, for all x belonging to the same isostable, Definition [2] 
imposes the same value |si(x) | in the decomposition of the flow (I12|) and the same asymptotic 
behavior ffT5j) . Note that without the multiplication by the increasing exponential e~ ait , one 
would have X T = B(x*) Vr since <f)(t, x) — x* — > as t — > oo for all x e £>(x*). 

Except for the case of multiple eigenvalues, for which vi might not be unique (see Remark 
[2]), the isostables are uniquely defined through Definition [2J Uniqueness is also related to 
the fact that the Koopman operator has a unique eigenfunction si (x) which is continuously 
differentiable in the neighborhood of the fixed point. 

Remark 1 (Unstable fixed point). Definition [2] is easily extended to unstable fixed points 
characterized by Oj > o\ > for all j. Indeed, the isostables are still given by Definition 
El where the limit t — > oo is replaced by t — > — oo, that is, one considers the flow 0(— t, x) 
induced by the (stable) backward-time system. In this case, the isostables are related to the 
unstable eigenfunction si(x) of the Koopman operator. 

Remark 2 (Multiple eigenvalues). When the eigenvalue Ai has a multiplicity m > 1, the 
fixed point is either a star node (m linearly independent eigenvectors) or a degenerate node 
(m linearly dependent eigenvectors). In the case of a star node, Definition [2J is not unique 
since it depends on the direction of the eigenvector V! (in other words, a C 1 eigenfunction 
of the Koopman operator corresponding to the eigenvalue Ai is not unique). Actually, v x 
should be replaced in Definition [2] by any linear combination of m orthonormal eigenvectors 
of Ai, a situation where the isostables lying in the vicinity of the fixed point correspond to 
cylindrical hypersurfaces whose intersection with the hyperplane spanned by the eigenvectors 
of Ai is a hypersphere. In the case of a degenerate node, the asymptotic evolution toward the 
fixed point is dominated by the (slowest) term si(x)v! t m_1 exp^t). Then, the increasing 
exponential exp(— ait) in Definition [2J must be replaced by t l ~ m exp(— <J\t). 



Some remarks on the isostables 



Equivalent definitions for excitable systems. In [24[], the authors considered two- 
dimensional excitable systems characterized by a transient attractor (i.e. slow manifold) 
which attracts all the trajectories as they approach the fixed point. They defined the isosta- 
bles (they actually used the term "isochrons", see Section [V AI) as the sets of points that 
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converge to the same trajectory on the transient attractor. This definition is equivalent to 
Definition [2] since both impose that trajectories on the same isostable have the same asymp- 



totic behavior (see also Section TlV A[) . However, the definition of 24j is qualitative since no 
trajectory effectively reaches the transient attractor (which may even lose its normal stabil- 
ity property near a fixed point with complex eigenvalues). Also, it is valid only if the system 
admits a transient attractor induced by the slow-fast dynamics. In contrast, Definition [2] is 
more general and does not rely on the existence of a transient attractor. 



For systems with a slow (or center) manifold, the "projection manifolds" studied in [4|.|26| 
are related to the isostables. They are the sets of initial conditions for which the trajectories 
share the same long-term behavior on the slow manifold. In addition, they can be obtained 



through the normal form of the dynamics \&\. If the slow manifold is one-dimensional and 
if Ai is real, the projection manifolds are identical to the isostables. Otherwise, they do not 
exactly correspond to the isostables since they are not related to the slowest direction Vi 
only and are not of codimension-1. 

Isostables and flow. The flow 4>(At, ■) maps the isostable X T to the isostable X r+ At, for 
all At G R (as explained in the beginning of Section [TTJ) . Indeed, if x G T T , Definition [2] 
implies that 





lim e~ ait 

t— >oo 



<f>(t, x) - x* - 3? {vi e l(wi * +e) } e ai{t+T) 



for some 6 G O. Using the substitution t — t' + At, we have 

lim e- ait ' (t', <f>(At, x)) - x* - & ( Vl e t ^ t ' +e ' ) ) e ai{t ' +T+At ^ 

t'-»oo I. > 







with 9' = 9 + wiAt G 6, so that <f)(At, x) G 1 T +At- 

Local geometry near the fixed point. The isostables close to the fixed point have a ge- 
ometry similar to the isostables of the linearized dynamics, i.e. hyperplanes (for Ai G R) or 
cylindrical hypersurfaces (for Xi ^ R) (see Section Til Aj) . This follows from the fact that, in 
the vicinity of the fixed point, the flow (j!2j) and the flow induced by the linearized dynamics 
are (approximately) equal, so that their eigenfunctions s\ (x) have (approximately) the same 
value for ||x — x*|| <C 1 (see also (IA4I) in Appendix lAl). 

Invariant fibration. When the eigenvalue Ai is real, the isostables are the invariant fibers 
of the 1-dimensional invariant manifold V defined as the trajectory associated with the slow 
direction Vx (i.e. the transient attractor in the case of slow-fast systems). Given their 
local geometry, it is clear that the isostables near the fixed point are the fibers defined by 
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the splitting N (BTV, where N = span{v 2 , . . . , v n } and TV = span{vi}. Moreover, it 
follows from the invariance property of the isostables that this local fibration is naturally 
extended to the whole invariant manifold V by backward integration of the flow. Provided 
that 02 < 0i, the normal hyperbolicity of V implies that the isostables are characterized by 
smoothness properties and persist under a small perturbation of the vector field j^, [l^ . In 
addition, this description also implies the uniqueness of the concept of isostables. Note that 
Definition [2] is recovered in |{|, Theorem 3, and corresponds to the property that the points 
on the same fiber converge to a trajectory on V with the fastest rate. 

When Ai is complex, however, the isostables cannot be interpreted as the fibers of an 
invariant manifold. In this case, their only definition is in terms of an eigenfunction of the 
Koopman operator. 

Extension to other eigenf unctions. The isostables X T are related to the first eigenfunc- 
tion Si(x) of the Koopman operator, but the concept can be directly generalized to other 
eigenfunctions. Namely, the sets I {]), j G J = {1, . . . , n} \ {i G N|A.j_i ^ M}, are obtained 
by considering the level sets of |s 3 -(x)|. The extension is useful to derive an action-angle 
coordinates representation of the system, to perform a global linearization of the dynamics 
(see Section l"VBj) . or to compute the (un)stable manifold of an attractor. 

The intersection between the sets I^, with j < j, is defined as the generalization of 
Definition [2] 



n x% = \ x g s« 

jej 



36j G Bj s.t. 



1 - —o—t 

hm e 3 

t— >oo 



0(t,x) -x* - J2^{-v j e i ^ t+e ^}e^ t+ri ^ 



(14) 



with Bj = {0, ?r} if Xj G R and 6,- = [0, 2tt) if Xj (£ E. When r (j) = oo for all j < j G J, 
f lT4|) is equivalent to Definition [21 so that it can be interpreted as an isostable for the system 
restricted to the invariant manifold M- = n, c -7 j,;2^1 . (The manifold M-- is associated 
with the fast directions Vj, j = j , . . . ,n.) In addition, if Aj G M, 0141) defines a codimension- 
j invariant fibration of the invariant manifold Vj = f]jej j>J-^vO)=oo' (The manifold Vj is 
associated with the slow directions v,-, j = If Vj is a slow manifold, then the 

fibration f[T4|) corresponds to the projection manifolds considered in (4, 26]. Note that the 



family of manifolds Vj generalizes the notion of slow manifold observed for systems with 

14 



slow-fast dynamics. 



III. LAPLACE AVERAGES 



In this section, we show that the isostables can be obtained through the computation of 
the so-called Laplace averages. The Laplace averages of a scalar observable / : MJ 1 H- C are 
given by 

/ A *(x) = lim - A/°0t)(x)e- At dt, (15) 

T-s-co 1 Jo 

with </>t(x) = (j)(t, x) and A G C. (The observable / has to satisfy some conditions which 
ensure that the averages exist.) When it exists and is nonzero for some A and /, the Laplace 
average /^( x ) corresponds to the eigenfunction of the Koopman operator associated with 
the eigenvalue A [20|. Indeed, one easily verifies that 



U*ft(x) = lim i f(/o^)(x)e" A! di 

f 

(/o0 i )(x)e- A *dt 



At' v 1 rT+t 

e xt hm — 



T^oo T Jt 



^ A*(x) 



where the second equality is obtained by substitution. For systems with a stable fixed point, 
the Laplace average (x) corresponds (up to a scalar factor) to the eigenfunction si(x), 



and is therefore related to the concept o 



extension of the Fourier averages 
limit cycles. 
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isostable. In addition, the Laplace averages are an 



2l| that were used in 18j to compute the isochrons of 



Remark 3. Instead of fll5p . the generalized Laplace averages 20| 



^ W = t 1 ^ T Jo ( (/ ° ^ (X) - /(X * } - £ & 



x)e Afc * e~ Xjt dt 



must be considered to obtain other eigenfunctions Sj (x), j > 2, and the associated sets IP' 



considered in ( I14|) . However, their computation is delicate since it requires a very accurate 
computation of the other (generalized) Laplace averages /? (x), k < j, and goes beyond the 
scope of the present paper. 
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A. The main result 



The exact connection between the Laplace averages and the isostables is given in the 
following proposition. 

Proposition 1. Consider an observable f G C 1 such that /(x*) = and (V/(x*), vx) ^ 0. 
Then, a unique level set of the Laplace average |/f | corresponds to a unique isostable. That 
zs > I/ai( x )I = I/ai( x ')I; with x G T T and x' G T T i , if and only if r = r' . In addition, 



T — T 



log 



KM) 



Proof. If x belongs to the isostable X T , one has, for some 9 G 0, 



lim e~ ait 

t— >oo 



(/ o t )(x) - /(x*) - (v/(x*), 3? { Vl e*^)}) e^ i+T 

'(V/(x*),&(x) - x*> + o(||&(x) - x*||) - (v/(x*),^{v ie ^ +e )}\^^) 
< ||V/(x*)|| lim e~ ait h t (x) - fx* + ^| Vl e i(wi * +e) } e CTl( ' +r) 



lim e~ CT1 * 

t— >oo 



+ lim e" CT1 *o 

t— »oo 



E 



^(x)---4»(x) v fcl ... fe „ e (*^i+-+*»^)* 



{fei,-,fcn}eNg- 

felH hfcn>l 







(16) 



The first equality is obtained through a first-order Taylor approximation, the inequality 
results from the Cauchy-Schwarz inequality and the expression of the flow (TT21) . and the last 
equality is implied by Definition [2J Then, it follows from f)16p that 

1 rT 1 rT 

'->oo ' 

1 r r 



lim - / (/o0 t )(x)e- Alt dt- lim - / (f(x*) + (V/(x*), 3? { Vl e^ f +^}) e CT ^>) e"* 1 ' <ft 
(/ o 0,)(x) - /(x*) - (v/(x*), 3? {vi e*^*^)}) e-^ 



< lim — / e 

- T^oo T Jo 

or equivalently, given (Il5j) and since /(x*) = 
1 r T 



dt = 



/X(x) =1^1 (/(x*) + (v/(x*),3ft{v ie ^ +e )})e-^) e- Al 'dt 
= lim I / T (v/(a 

T^oo 1 Jo \ 



= lim — ( / T (V/(x*), Vl ) e <T1T+ie dt+ / T (V/(x*), vj) e^^^t+e) df \ ^ 
t->oo 2T yjo jo / 

If Ai G K, one has Wi = 0, Vi = vf, and e ie = e~ ie (since 6> G = {0, 7r}). Then, it follows 

from (jnp that 

1 /" T 



/^(x) = lim - / (V/(x*), Vl > e n *»dt = (V/(x*), Vl > e 



(18) 



1(3 



and 

\f* Xi (x.)\ = \(Vf(x*),v 1 )\e^, AiGR. (19) 
If Ai ^ K, wi 7^ implies that the second term of (fTTl) is equal to zero, which yields 

l^(x)|= l(V/(X 2 * ),Vl) ' e^, A^R. (20) 

For x' G Xj-i, the inequalities (TT9"j) or (I2"01 still hold (with r replaced by r'), so that the result 
follows provided that (V/(x*),Vi) ^ 0. □ 

Remark 4 (Unstable fixed point and multiple eigenvalues (see also Remarks [TJ and [2])). (i) 
For unstable fixed points with <jj > o\ > for all j, the isostables are the level sets of the 
Laplace averages |/*^J computed for backward-in-time trajectories 4>{—t, •). 

(ii) In the case of a star node (e.g. with a real eigenvalue of multiplicity m), the isostables 
obtained through the Laplace averages depend on the choice of the observable /, which 
can have a nonzero projection (V/(x*), v^), j — 1,. . . ,m, on several eigenf unctions of the 
Koopman operator associated with the eigenvalue Ai. However, a unique family of isostables 
is obtained by considering the level sets of \jYlk n =i(f\ 1 ,k)' 2 > where fl i k denotes the Laplace 
average for an observable fk that satisfies (V/fc(x*), Vj) = for all j G {1, . . . , m} \ {k}. 

(iii) In the case of a degenerate fixed point (eigenvalue of multiplicity m), the isostables 
are computed with the Laplace averages, but the exponential exp(— Ait) in fTTol) must be 
replaced by t x ~ m exp(— Ait). 



B. Numerical computation of the Laplace averages 

Proposition [T] shows the strong connection between the isostables and the Laplace aver- 
ages, a result which provides a st raig htforward method for computing the isostables. Simi- 



larly to the method developed in [18| , the computation of isostables is realized in two steps: 
(i) the Laplace averages are computed (over a finite time horizon) for a set of sample points 
(distributed on a regular grid or randomly); (ii) the level sets of the Laplace averages (i.e. 
the isostables) are obtained using interpolation techniques. The proposed method is flexible 
and well-suited to the use of adaptive grids, for instance. In addition, the averages can be 
computed either in the whole basin of attraction of the fixed point or only in regions of 
interest. 
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It is important to note that the computation of the Laplace averages involves the multi- 
plication of the very small quantity (/ o <^)(x) with the very large quantity exp(— A]i), as 
t — > oo. When the trajectory approaches the fixed point, the relative error of the integration 
method implies that the (numerically computed) quantity (/ o 0t)(x) does not compensate 
exactly the value exp(— Ait), and the computation becomes numerically unstable. Therefore, 
a high accuracy of the numerical integration scheme and a reasonably small time horizon T 
are required for the computation of the Laplace averages. 

In spite of the numerical issue mentioned above, an algorithm based on a straightforward 
calculation of the Laplace averages produces good results. However, it is improved if one can 
avoid computing the integral. Toward this end, we remark that evaluating the integral (|15p 
is not necessary when Ai is real, since the integrand converges to a constant value. When Ai 
is complex, we consider the successive iterations of the discrete time-Xi map 4>(Ti, •), with 
T\ = 2ir/oj\. The result is summarized as follows. 

Proposition 2. (i) Real eigenvalue \\. Consider an observable f G C 1 that satisfies /(x*) = 
0. Then, the Laplace average (x) corresponds to the limit 

ft(x) = lime-^(/o^)(x). 

(ii) Complex eigenvalue \\. Consider two observables fx G C 1 and f 2 G C 1 that satisfy 

/i(x*) = / 2 (x*) = 
|(V/i(x*),a)| = |(V/ 2 (x*),b)|^0 
(V/i(x*),b) = (V/ 2 (x*),a) = 

with a = 3?{vi} and h = — S{vi}. Then the Laplace average |/JJ (x)| of an observable 
f G C 1 is proportional to the limit 

| f* Xi (x) | oc n lim v / ((Ao0 nTl )(x)) 2 + ((/ 2 o0 nTl )(x)) 2 , 

neN 

with T\ = 2-k/ui. 

Proof, (i) Real eigenvalue \\. Since /(x*) = 0, the result follows from f)16p and (I18p . 
(ii) Complex eigenvalue X\. Provided that /(x*) = 0, (TT5|) implies that 

lirn e- ffinTl (/ofe)(x) = (V/(x*), 3? Wie 10 }) e^ T = (V/(x*), a cos(9) + b sin(0)> e^ T 
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and, in particular, 

Jim e- CT1 " Tl (/io0 nTl )(x) = cos(£)(VA(x*),a>e CT1T , 
Jim e- CTl " Tl (/ 2 o0 nri )(x) = sin(0)(V/ 2 (x*),b)e^. 

Then, one has 

Jim e-^VCCA ° ^)(x)) 2 + ((/ 2 o nTl )(x)) 2 = | (V/xM.a) |e CT - 

and it follows from (I2U|) that the limit is proportional to \f^ i (x*)\ — with the factor of pro- 
portionality 2| (V/i(x*), a) / (V/(x*), vi) |. □ 

Proposition [2] implies that the isostables can be computed as the level sets of particular 
limits. While the straightforward computation of the Laplace averages (1151) is characterized 
by a rate of convergence T _1 , the computation of these limits is characterized by an expo- 
nential rate of convergence. Hence, the results of Proposition [2] are of great interest from a 
numerical point of view, and it is particularly so since the numerical instability imposes an 
upper bound on the finite time horizon T. 

IV. APPLICATIONS 

The concept of isostables of fixed points is now illustrated with some examples. These 
examples show that the framework is coherent and general, coherent with the equivalent 
definition of isostable for excitable systems and general since it is not limited to the particular 
class of excitable systems. 

A. The excitable FitzHugh-Nagumo model 

The concept of isostables is primarily motivated by the reduction of excitable systems 
characterized by slow-fast dynamics. In this case, the points on the same isostable X r share 
the same asymptotic behavior on the stable slow manifold. 

In this context, we compute the isostables for the well-known FitzHugh-Nagumo model 

QQ 

v = —w — v(v — l)(v — a) + I , 
w = e(v — •yw) , 
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which admits an excitable regime with a stable fixed point (x* = (v*,w*), with v* = w*) 
for the parameters / = 0.05, e = 0.08, 7 = 1, and a = {0.1,1}. The eigenvalues (of the 
Jacobian matrix at the fixed point) are either real (e.g., a = 1) or complex (e.g., a = 0.1). 
We consider both cases in the sequel. 



In 



24( , the isostables were computed for the FitzHugh-Nagumo model through the back- 



ward integration of trajectories starting in a close neighborhood of the stable slow manifold 
(or transient attractor). Here, we obtain the same results using a forward integration method 
based on the computation of the Laplace averages. 



1. Real eigenvalues (a = 1) 

The Laplace averages are computed according to the result of Proposition r2]d), with the 
observable f(v, w) = (v — v*) + (w — w*). The level sets (isostables) are represented in Figure 




Figure 4: The level sets of the Laplace averages |/f | are the isostables (black curves) of the fixed 
point (red dot). In the neighborhood of the fixed point, the isostables are parallel to the direction 
v 2 ~ (—1,0.11) (red arrow). Two trajectories with an initial condition on the same isostable 
((—0.0303, —0.5152) for the solid curve, (1.7879, —0.8182) for the dashed curve) synchronously 
reach the same isostable after a time r — r' ~ 12. They also reach the stable slow manifold 
(transient attractor) (green curve) synchronously. (The averages are computed on a regular grid 
100 x 100, with a finite time horizon T = 50; the black dotted-dashed curves are the nullclines.) 



One first verifies that the isostables are parallel to the eigenvector V2 in the neighborhood 
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of the fixed point. In addition, two trajectories with an initial condition on the same isostable 
synchronously converge to the fixed point. For instance, two trajectories that start from the 
same level set |si(x')| = 1.74 synchronously reach the level set |si(x)| = 0.17 after a time 
t — t' ~ 12. This observation confirms the result of Proposition 1, since 



— log 

0"! 



Si(Xl 



0.17 
log —— & 12 . 



-0.1933 °1.74 



The system admits an unstable slow manifold (transient repeller), which corresponds to 
a stable slow manifold (transient attractor) for the backward-time system. The unstable 
slow manifold lies in the highly sensitive region v < 0, w ~ —0.3 characterized by a high 
concentration of isostables. If a trajectory lying near the fixed point on an isostable X T is 
perturbed beyond the unstable slow manifold, it will reach an isostable Z T ', with r' <C r. As 
a consequence, the trajectory will not immediately converge toward its initial position near 
the fixed point but will exhibit a large excursion in the state space, whose duration is given 
by t - r' > 1. This is a characteristic phenomenon of slow-fast excitable systems. Note 
that for slow-fast asymptotically periodic systems, a high concentration of isochrons is also 
observed near the unstable slow manifold 



2. Complex eigenvalues (a = 0.1J 

The Laplace averages are computed according to the result of Proposition [2](ii), with the 
observables fi(v,w) = b%{v — v*) — b\(w — w*) and f 2 (v,w) = a 2 (v — v*) — ai(w — w*), 
a = (01,02), b = (61,62)- The level sets (isostables) are represented in Figure EJ We 
verify that the isostables are ellipses in the neighborhood of the fixed point (Figure 0(b)). 
In addition, two trajectories with an initial condition on the same isostable synchronously 
converge to the fixed point (Figure 0(a)). For instance, two trajectories that start from the 
same level set |si(x')| = 0.10 synchronously reach the level set |si(x)| = 0.051 after a time 
t — t' ~ 16. This observation confirms the result of Proposition 1, since 



— log 



Si X 



Sl X' 



1 °- 051 ~ir 
log ~ lb 



-0.041 0.10 



As in the case Ai real, the system admits a unstable slow manifold (region v < and 
w ~ 0) characterized by a high concentration of isostables. 
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(a) (b) 



Figure 5: The level sets of the Laplace averages |/^ | are the isostables (black curves) of the 
fixed point (red dot), (a) Two trajectories with an initial condition on the same isostable 
((0.7688,-0.5779) for the solid curve, (—0.1960,-0.1558) for the dashed curve) synchronously 
reach the same isostable after a time r — r' ~ 16. (The averages are computed on a regular 
grid 100 x 100, with a finite time horizon T = 250, that is, with 11 iterations of the time-Ti 
map; the black dotted-dashed curves are the nullclines.) (b) In the neighborhood of the fixed 
point, the isostables are ellipses. The arrows represent the vectors a = 5R{vi} ~ (0.96,0.03) and 
b = — 9{vi} ps (0, 0.27). (The averages are computed on a regular grid 50 x 50). 

B. The Lorenz model 

The framework developed in this paper is not limited to two-dimensional excitable mod- 
els, but can also be applied to higher-dimensional models, including those which are not 
characterized by slow-fast dynamics. For instance, we compute in this section the isostables 
of the Lorenz model 

x'i = a(x 2 - xi) , 

X2 = %i{p - x 3 ) - x 2 , 

x 3 = xix 2 - bx 3 . 

With the parameters a — 10, p — 0.5, 6 = 8/3, the origin is a stable fixed point with a real 
eigenvalue Ai. Several isostables are depicted in Figure El They are the two-dimensional 
level sets — i.e., the isosurfaces — of the Laplace averages fl computed for the observable 
f(xi,X2,Xs) = X\ + £2 + 23. Note that the isostables are approximated by a plane in the 
vicinity of the fixed point. 

When the parameter p exceeds the critical value p — 1, the origin becomes unstable and 
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7 "5 



Figure 6: The isostables can be computed for three-dimensional models, including those which 
are not characterized by slow-fast dynamics (in this case, the Lorenz model). Four isostables are 
represented, which are the level sets of the Laplace averages |/^J G {0.5, 1, 1.5,2}. (The averages 
are computed on a regular grid 75 x 75 x 75, with a finite time horizon T = 20; the red dot 
corresponds to the fixed point.) 

two stable fixed points (±x*, ixJj, Xg) appear. Since these fixed points are characterized 
by the same eigenvalues, their isostables can be obtained through the computation of a 
single Laplace average . In Figured the isostables are computed for the value p = 2, a 
situation characterized by a complex eigenvalue Ai. Note that the isostables are cylinders 
in the vicinity of the fixed point. In addition, the level set |/^J — > oo corresponds to the 
separatrix between the two basins of attraction (i.e. the stable manifold of the fixed point 
at the origin). 

V. DISCUSSION 

In this section, we discuss some topics related to the concept of isostables. Through 
the Koopman operator framework, we claim that the notion of isostables is different from 
but complementary to the known notion of isochrons. Isostables and isochrons define a 
set of action-angle coordinates and are related to a global linearization of the dynamics. 
In addition, we briefly show that the isostables are the level sets of a particular Lyapunov 
function for the fixed point dynamics. 



23 



-5 4 



2 X20 



-2 -4 



Figure 7: The level sets of the Laplace averages \ f% | € {1, 2, 3, 4, 5} represent five isostables of the 
two stable fixed points. (The averages are computed on a regular grid 50 x 50 x 50, with a finite 
time horizon T = 15; the red dot corresponds to the (visible) stable fixed point.) 



A. Isostables vs. isochrons 



The isostables are the sets of points that approach the same trajectory when they con- 
verge toward the fixed point. Similarly, in the case of asymptotically periodic systems, the 
isochrons are the set of points that converge toward the same trajectory on the limit cycle 



321 ] . It follows that isostables (of fixed points) and isochrons (of limit cycles) are conceptu- 
ally related. However, these two concepts are also characterized by intrinsic differences and 
turn out to be complementary. 

The difference between isostables and isochrons can be understood through the framework 
of the Koopman operator. The isostables have been defined as the level sets of the absolute 



value of the Koopman eigen: 
cycles were computed in 



'unction |si(x)| (Section III Bj) . In contrast, the isochrons of limit 
by using the argument of a Koopman eigenfunction. Similarly, 
the isochrons of fixed points (characterized by a complex eigenvalue Ai) can be defined as 
the levels sets of the argument Zsi(x). These sets (also called isochronous sections) are 
well-known and usually defined as the sets invariant under a particular return map (i.e. the 
discrete map 0(T 1; •) considered in Proposition [2]). Also, their existence, which is not trivial 
in the case of weak foci (i.e. imaginary eigenvalues) or nonsmooth vector fields, has been 
investigated in 0, In the case of linear systems, the isochrons correspond to radial 
lines that intersect at the fixed point (see Figure E(b)). For nonlinear systems, they are 
characterized by a more complex geometry, even though they are still radial lines in the 
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vicinity of the fixed point (see Figure [S]). 

Isostables and isochrons appear to be two different but complementary notions. On 
one hand, the isostables are related to the stability property of the system and provide 
information on how fast the trajectories converge toward the attractor. On the other hand, 
the isochrons are related to a notion of phase and provide information on the asymptotic 
behavior of the trajectories on the attractor. Given (II ip . the isostables are related to the 
property 

±\ ai (M*))\ = °iMM*))\ (2i) 

while the isochrons are characterized by 

^Zs 1 (0 4 (x))=u; 1 . (22) 

In the case of fixed points, it is clear that the isochrons are not relevant to characterize the 
synchronous convergence of the trajectories, a fact that stresses the importance of considering 
the isostables instead. 

0.2 
0.1 

W 

-0.1 
-0.2 

-0.2 -0.1 0.1 0.2 

Figure 8: For a fixed point with a complex eigenvalue Ai, the isostables (black curves) and the 
isochrons (red curves) of the fixed point are the level sets of |si(x)| and Zsi(x), respectively. 
In the vicinity of the fixed point, the isostables are ellipses and the isochrons are straight lines. 
(The numerical computations are performed for the FitzHugh-Nagumo model, with the parameters 
considered in Section [IV A 2t the blue dot represents the fixed point.) 
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B. Action-angle coordinates and global linearization 



For a two-dimensional dynamical system which admits a spiral sink (two complex eigen- 
values), the families of isostables and isochrons provide an action-angle coordinates repre- 
sentation of the dynamics. More precisely, fl2~Tj) and fl2~2"j) imply that, with the variables 
(r, 9) = (|si(x)|, Zsi(x)), the system is characterized by the (action-angle) dynamics 



r = a\r 

9 = 0)! 

in the basin of attraction of the fixed point. For systems of higher dimension, the action- angle 
dynamics are obtained with several Koopman eigenfunctions, i.e. (rj, 9j) = (|sj(x)|, Zsj(x)) 
leads to fj = crjTj, 9j = ujj. Note that this was also shown in Section Hi A 21 in the case of 
linear systems with a spiral sink. 

When expressed in the action-angle coordinates, the dynamics become linear. This is 



in agreement with the recent work 



151 ] showing that a coordinate system which linearizes 



the dynamics is naturally provided by the eigenfunctions of the Koopman operator (see also 
Appendix [A]). Namely, in the new variables yj = Sj(x), the system dynamics are given by 

fx, 



d_ 

dt 



\y„J 

Moreover, the linear change of coordinates 



\0 A 2 J 



\ynj 



v 



\Z n ) 



\VnJ 



(23) 



where the columns of V are the eigenvectors Vj of the Jacobian matrix J at the fixed point, 
leads to the linear dynamics 



d 
dt 



\ z «/ 



For the two-dimensional FitzHugh-Nagumo model, the coordinates (z%, z^) are represented in 
FigureEland are equivalent to the action-angle coordinates (r, 9) (FigureEJ). They correspond 
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to Cartesian coordinates in the vicinity of the fixed point, where the linearized dynamics are 
a good approximation of the nonlinear dynamics (see also (1A3|) in Appendix . But owing 
to the nonlinearity, the coordinates are deformed as their distance from the fixed point 
increases. The comparison between these coordinates and regular Cartesian coordinates 
therefore appears as a measure of the system nonlinearity. 

In the case of two-dimensional systems with a stable spiral sink, the derivation of action- 
angle coordinates and the global linearization are obtained through the isostables and the 
isochrons, that is, with only the first Koopman eigenfunction Si(x). For more general 
systems, global linearization involves several Koopman eigenf unctions Sj(x) (see 15J for a 
detailed study), which can be obtained through the generalized Laplace averages (see Remark 
E]). In the context of model reduction, or when the dynamics are significantly slow in one 
particular direction, the first eigenfunction — related to the isostable — is however sufficient 
to retain the main information on the system behavior. 




Figure 9: The coordinates z\ (black curves) and Z2 (red curves) correspond to Cartesian coordinates 
in the vicinity of the fixed point but are deformed when far from the fixed point. (The numerical 
computations are performed for the FitzHugh-Nagumo model, with the parameters considered in 
Section HV A 21 the blue dot represents the fixed point.) 



C. Lyapunov function and contracting metric 

As a consequence of the linearization properties illustrated in the previous section, the 
Koopman eigenfunctions — and in particular the isostables — can be used to derive Lyapunov 
functions and contracting metrics for the system. 
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In the particular case of two-dimensional systems with a spiral sink, the isostables are the 
level sets of the particular Lyapunov function V(x) = |si(x)| (see Figure [TU1 for the FitzHugh- 
Nagumo model). Indeed, (J2H) implies that V(x) = 0iV(x) < Vx G 6(x*) \ {x*} and one 
verifies that V(x*) = 0. This function is a special Lyapunov function of the system, in the 
sense that its decay rate is constant everywhere. (Note that the function V = log(|si(x)|)/a"i 
satisfies V = — 1 but with V(x*) = — oo.) 




Figure 10: The function V = |si(x)| is a particular Lyapunov function for the system (here, the 
FitzHugh-Nagumo model with the parameters considered in Section IIV A2|) . One verifies that 
the function decreases with a constant rate along a trajectory (black curve). Note also that the 
unstable slow manifold (region v < 0, w ~ 0) is characterized by a line of maxima of the Lyapunov 
function. 

In addition, the isostables are related to a metric which is contracting in the basin of 
attraction of the fixed point. Namely, the distance 

d(x,x') = |si(x) - si(x')| 

is well-defined and ( II 1|) implies that 

^d(&(x),&(x')) = M(x,x') < 0, Vx + x' G B(x*) . 

For more general systems that admit a stable fixed point, the function V(x) = |si(x)| is still 
decreasing along the trajectories, but V(x) = does not imply x = x* (V is zero on the 
whole isostable X r=00 that contains the fixed point). However, the function can be used with 
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the LaSalle invariance principle. To obtain a good Lyapunov function, several Koopman 
eigenfunctions must be considered. For instance, the function 

i/p 

v(x)=lj>j 

with the integer p > 1, satisfies 



v(x) = £ l*i(x)l p E^M*)I P < o-iV(x) 

and V(x) = iff x = x*. In addition, a contracting metric is given by 

( i/p 

d(x,x') = ( E |si(x) - Sj(x!) 
and one has 

d -d(0 t (x),^(x')) < M^x') , Vx,x' G B(x*) . 



It follows from the above observations that the existence of stable eigenfunctions of the 
Koopman operator is closely related to the stability properties of the attractor. Therefore, 
the Koopman operator framework could potentially lead to a new method for the global 
stability analysis of nonlinear systems. 



VI. CONCLUSION 



In this paper, the well-known phase reduction of asymptotically periodic systems has 
been extended to the class of systems which admit a stable fixed point. In the context 
of the Koopman operator framework, the approach is not restricted to excitable systems 
with slow-fast dynamics but is valid in more general situations. The isostables required for 
the reduction of the dynamics, which correspond is some cases to the fibers of a particular 
invariant manifold of the system, are interpreted as the level sets of an eigenfunction of the 
Koopman operator. In addition, they are shown to be different from the concept of isochrons 
that prevails for asymptotically periodic systems. Beyond its theoretical implications, the 
framework also yields an efficient (forward integration) method for computing the isostables, 
which is based on the estimation of Laplace averages along the trajectories. 

The reduction of the dynamics through the Koopman operator framework leads to an 
action-angle coordinates representation that is intimately related to a global linearization 
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of the system. More precisely, the proposed reduction procedure is nothing but a global 
linearization of the system where only one direction of interest is considered, which retains 
the main information on the system behavior (i.e. the slowest direction). In this context, 
the isostables — related to the action — or the isochrons — related to the angle — used for the 
reduction are particular objects involved in the global linearization process. Given this 
relation between reduction methods and linearization, research perspectives are twofold. 
On the one hand, convenient (generalized) Laplace average methods could be developed for 



linearization purposes (e.g. computation of the isostables of limit cycles but for the 

computation of (un)stable manifolds as well. On the other hand, the Koopman operator 
framework can be further exploited for the reduction of more general dynamical systems 
(e.g. chaotic systems). 
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Appendix A: Spectral decomposition of the Koopman operator 

In this appendix, we derive the expansion fllOp of an observable onto the eigenf unctions of 
the Koopman operator. Consider the change of variable s : x y y, with jjj = Sj-(x), where 
Sj is an eigenfunction of the Koopman operator. It follows that s(x*) = and, given (1111) . 
the dynamics is linearized in the y variable, i.e. yj = XjVj- According to the linearization 
Poincare theorem 8], the transformation s is analytic since the vector field F is analytic and 
the eigenvalues are nonresonant (and provided there is no unstable fixed point in B(x*)). If 
an observable / is analytic, the Taylor expansion of /(s _1 (y)) around the origin yields 

f(s-\y)) = /(x*) + Vf(x^) J s - iy + V J^HJ.-* y + V ± * 



k 

x* 



H sr iy + h.o.t. , 
(Al) 

where J s -i is the Jacobian matrix of s -1 at the origin (i.e. Js- 1 ,^ — ds^ 1 /dyj(0)), H is 
the Hessian matrix of / at x* (i.e. Hij = d 2 f/(dxidxj)(x.*)), and H s i is the Hessian 
matrix of s^ 1 at the origin (i.e. H s -i ^ = d 2 s^ 1 / (dyidyj)(0)). Using the relationship y = 
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(si(x), . . . , s„(x)), we can turn the expansion (lAlj) into an expansion of / onto the products 
of the eigenfunctions Sj. For a vector- valued observable f , we obtain 



ffxl 



E V fel-fe» S l 



fcl , 



(A2) 



{fel,...,fc n }€N« 



with the (first) Koopman modes 



Vfe 1 ...fc n — < 



f(x*) 



dsZ 1 



dy 3 



yy 

t^ti dx k d Xl 

1 EE 



dx k dxi 



o % 

dsr 1 



E 



0f 



fc=l ® X k 



o % 



1 " df 



kj = l, h = Vi 7^ j , 

ki = kj = 1 , k r = Vr 7^ {i, j} 

jfci = 2, fcj- = OVj ^ i. 



fc=i 

The other (higher-order) Koopman modes can be derived similarly from flAll) . Since the 
eigenfunctions satisfy (I 111) , the relationship fllOp directly follows from (1A2j) . 
For the observable /(x) = x, the Koopman modes are given by 



v fcl ... fc 

n 



fci!...A; n ! d k ^yi - ■ -d k ™y n ' 

In particular, the eigenvectors of the Jacobian matrix J of F (i.e. Vj = v^...^, with fcj = 1, 
ki = Wi j) correspond to 

9s ~ 



7 % 

and one has J s -i = V, where the columns of V are the eigenvectors Vj. It follows that the 
variables z introduced in (12"3"1) satisfy z = J s -iy so that flAll) implies 



x = x* + z + o z . 



In addition, the derivation of y = s(s *(y)) at the origin leads to 



(A3) 



= ( V Si (x*) 



ds~ 



V Sl (x*), Vj c 



Therefore, the gradient Vsj(x*) is the left eigenvector of J (associated with the eigenvalue 
A,) and one has 



*(x) = (x - x*, Vs?(x*)) + o(||x - x*||) = (x - x*, v<> + o(||x - X* 



(A4) 
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which implies that, for ||x — x*|| <C 1, the eigenfunction Sj(x) is well approximated by the 
eigenfunction of the linearized system. 
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